load("SimNets.RData")
load("networkClean.RData")
source('NetSims/totalPathReduction.R')
library(sna)
library(network)

nodes <- nodeTitles

amat <- symmetrize(amat)

grange <- range(which(substr(nodes,1,2)=="g_")) # regulations
rrange <- range(which(substr(nodes,1,2)=="r_"))	# articles

verts <- c(rrange[1]:rrange[2],grange[1]:grange[2])
ArtNet <- amat[verts,verts]

artIsol <- (rrange[1]:rrange[2])[which(degree(ArtNet)[1:length(rrange[1]:rrange[2])] == 0)]

nodes <- nodes[-artIsol]
amat <- amat[-artIsol,-artIsol]

doti <- which(nodes =="a_departmentoftheinterior")
doi <- which(nodes == "a_departmentofinterior")
amat[doi,] <- amat[doi,]+amat[doti,]
amat[,doi] <- amat[,doi]+amat[,doti]

nodes <- nodes[-doti]
amat <- amat[-doti,-doti]

#write.csv(cbind(nodeTitles,""),"nodeList.csv",row.names=F)

library(colorRamps)
require(sna)

grange <- range(which(substr(nodes,1,2)=="g_")) # regulations
jrange <- range(which(substr(nodes,1,2)=="j_"))	# journals
rrange <- range(which(substr(nodes,1,2)=="r_"))	# articles
srange <- range(which(substr(nodes,1,2)=="s_")) # scientists
frange <- range(which(substr(nodes,1,2)=="f_")) # funders
arange <- range(which(substr(nodes,1,2)=="a_")) # agencies
crange <- range(which(substr(nodes,1,2)=="c_")) # committees
urange <- range(which(substr(nodes,1,2)=="u_")) # affiliations

# Node names
gnames <- nodes[grange[1]:grange[2]] # regulations
jnames <- nodes[jrange[1]:jrange[2]]	# journals
rnames <- nodes[rrange[1]:rrange[2]]	# articles
snames <- nodes[srange[1]:srange[2]] # scientists
fnames <- nodes[frange[1]:frange[2]] # funders
anames <- nodes[arange[1]:arange[2]] # agencies
cnames <- nodes[crange[1]:crange[2]] # committees
unames <- nodes[urange[1]:urange[2]] # affiliations
  
## Find Positions for Nodes
verts <- c(rrange[1]:rrange[2],grange[1]:grange[2],arange[1]:arange[2])
ArtNet <- amat[verts,verts]
set.seed(5)
xy <- gplot(ArtNet)
pos <- xy[1:length(c(rrange[1]:rrange[2],grange[1]:grange[2])),]
apos <- xy[-(1:length(c(rrange[1]:rrange[2],grange[1]:grange[2]))),]

# Agency positions at centroids
anodes <- (1:nrow(ArtNet))[-(1:length(c(rrange[1]:rrange[2],grange[1]:grange[2])))]
apos <- matrix(0,length(anames),2)
for(i in 1:length(anodes)){
	nodesi <- which(ArtNet[anodes[i],]==1)
	xi <- pos[nodesi,1]
	yi <- pos[nodesi,2]
	apos[i,] <- c(mean(xi),mean(yi))
}


# Colors based on agency
set.seed(95)
acols <- sample(colors()[1:100],length(anames))
gseq <- grange[1]:grange[2]
aseq <- arange[1]:arange[2]
gcols <- numeric(length(gnames))
for(i in 1:length(gseq)){
	ri <- gseq[i]
	gcols[i] <- acols[which(amat[ri,aseq]==1)]
}

anetplot <- anames[which(!is.na(apos[,1]))]
apcols <- acols[which(!is.na(apos[,1]))]
apos <- apos[which(!is.na(apos[,1])),]

ashort <- c("DoAg","DoEn","HHS","DHS","HUD","DoI","DoJ","DoL","DoT","EPA","GSA","JBR","EPA&DoT")



pdf("RIAArticleNet.pdf",height=20,width=20)
set.seed(5)
library(sna)
# use mixed color for articles cited by multiple agencys dgcols 
### Remove Article Isolates
vcols <- c(rep("white",rrange[2]-rrange[1]+1),gcols)
verts <- c(rrange[1]:rrange[2],grange[1]:grange[2])
nside <- c(rep(3,length(rnames)),rep(25,length(gcols)))
size <- c(rep(1.2,length(rnames)),rep(2,length(gcols)))
ArtNet <- amat[verts,verts]
gplot(ArtNet,coord=pos,vertex.col=vcols,edge.col=rgb(15,15,15,alpha=75,maxColorValue=255),gmode="graph",vertex.sides=nside,vertex.cex=size)
library(berryFunctions)
for(i in 1:length(ashort)){
textField(apos[i,1],apos[i,2],ashort[i],col=apcols[i],cex=3,border=apcols[i])
}
dev.off()

## Centrality Analysis
type <- c(rep("art",rrange[2]-rrange[1]+1),rep("reg",grange[2]-grange[1]+1))
#betCent <- betweenness(ArtNet)
degCent <- degree(ArtNet,cmode="indegree")
bcArt <- pathMembership(ArtNet,which(type=="art"),which(type=="reg"))
dcArt <- degCent[which(type=="art")]


## Top 5 articles
rnames[order(-dcArt)][1:5]
dcArt[order(-dcArt)][1:5]

load("simBC.RData")

# Geodesic distance and betweenness
allDist <- geodist(amat)$gdist==2

# Article Funders
fundMat <- amat[rrange[1]:rrange[2],frange[1]:frange[2]]
colnames(fundMat) <- fnames

# Article University
univMat <- allDist[rrange[1]:rrange[2],urange[1]:urange[2]]
colnames(univMat) <- unames

affilDat <- read.csv("AffiliationData.csv",stringsAsFactors=F)

fundDat <- read.csv("FunderData.csv",stringsAsFactors=F)
isdept <- fundDat$department

# Identify University
isuniv <- affilDat$university

# Analysis of funders
isnatnl <- fundDat$government
#isfound <- numeric(length(fnames))
#isfound[grep("foundation",fnames)] <- 1

# Identify depts and agencies
#isdept <- numeric(length(fnames))
#isdept[grep("environmental protection agency",fnames)] <- 1
#isdept[grep("department of energy",fnames)] <- 1
#isdept[grep("national oceanic",fnames)] <- 1
#isdept[grep("department of agriculture",fnames)] <- 1
#isdept[grep("centers for disease control",fnames)] <- 1
#isdept[grep("fish and wildlife",fnames)] <- 1


## Analysis of Funders (proof of concept -- dept or agency funded)
deptArt <- fundMat%*%cbind(isdept)>0

mean_diffs_universityF <- numeric(length(simNets))
for(i in 1:length(simNets)){
	bci <- simBCArt[[i]]
	bcArti <- bci
	mean1 <- mean(bcArti[deptArt==1])
	mean0 <- mean(bcArti[deptArt==0])
	mean_diffs_universityF[i] <- mean1-mean0
	if(i/25 == round(i/25)) print(i)
	}

mean1 <- mean(bcArt[deptArt==1])
mean0 <- mean(bcArt[deptArt==0])
mean_diff_universityF <- mean1-mean0

pdf("BoxPlotDept.pdf",height=3.5,width=4,pointsize=11)
par(las=1,mar=c(2,4,1,1))
boxplot(mean_diffs_universityF,ylab="Mean Difference in Shortest Paths")
abline(h=mean_diff_universityF,lwd=1.5,col="red")
pgt <- mean(mean_diffs_universityF>mean_diff_universityF)
p = 2*min(c(1-pgt,pgt))
legend("topleft",legend=paste("p = ",p,sep=""),bty="n")
dev.off()

mean_diffs_DeptF <- mean_diffs_universityF-mean_diff_universityF

## Look at top 10 funders
fundDeg <- apply(fundMat,2,sum)
top50funders <- fnames[order(-fundDeg)[1:50]]

# Identify major, non-gov funder
isnongov <- (fundDat$government==0)*(fundDat$department==0)

ngovArt <- fundMat%*%cbind(isnongov)>0

mean_diffs_universityNGF <- numeric(length(simNets))
for(i in 1:length(simNets)){
	bci <- simBCArt[[i]]
	bcArti <- bci
	mean1 <- mean(bcArti[ngovArt==1])
	mean0 <- mean(bcArti[ngovArt==0])
	mean_diffs_universityNGF[i] <- mean1-mean0
	if(i/25 == round(i/25)) print(i)
	}

mean1 <- mean(bcArt[ngovArt==1])
mean0 <- mean(bcArt[ngovArt==0])
mean_diff_universityNGF <- mean1-mean0

mean(mean_diffs_universityNGF>mean_diff_universityNGF)
mean_diffs_NGF <- mean_diffs_universityNGF-mean_diff_universityNGF

## Analysis of Articles
univArt <- univMat%*%cbind(isuniv)>0

mean_diffs_university <- numeric(length(simNets))
for(i in 1:length(simNets)){
	bci <- simBCArt[[i]]
	bcArti <- bci
	mean1 <- mean(bcArti[univArt==1])
	mean0 <- mean(bcArti[univArt==0])
	mean_diffs_university[i] <- mean1-mean0
	if(i/25 == round(i/25)) print(i)
	}

mean1 <- mean(bcArt[univArt==1])
mean0 <- mean(bcArt[univArt==0])
mean_diff_university <- mean1-mean0

pdf("BoxPlotUniv.pdf",height=3.5,width=4,pointsize=11)
par(las=1,mar=c(2,4,1,1))
boxplot(mean_diffs_university,ylab="Mean Difference in Shortest Paths")
abline(h=mean_diff_university,lwd=1.5,col="red")
pgt <- mean(mean_diffs_university>mean_diff_university)
p = 2*min(c(1-pgt,pgt))
legend("topleft",legend=paste("p = ",p,sep=""),bty="n")
dev.off()

mean_diffs_UnivU <- mean_diffs_university-mean_diff_university


## Analysis of Articles

# Identify depts and agencies
isudept <- affilDat$usgov

udeptArt <- univMat%*%cbind(isudept)>0

mean_diffs_university <- numeric(length(simNets))
for(i in 1:length(simNets)){
	bci <- simBCArt[[i]]
	bcArti <- bci
	mean1 <- mean(bcArti[udeptArt==1])
	mean0 <- mean(bcArti[udeptArt==0])
	mean_diffs_university[i] <- mean1-mean0
	if(i/25 == round(i/25)) print(i)
	}

mean1 <- mean(bcArt[udeptArt==1])
mean0 <- mean(bcArt[udeptArt==0])
mean_diff_university <- mean1-mean0

pdf("BoxPlotUDept.pdf",height=3.5,width=4,pointsize=11)
par(las=1,mar=c(2,4,1,1))
boxplot(mean_diffs_university,ylab="Mean Difference in Shortest Paths")
abline(h=mean_diff_university,lwd=1.5,col="red")
pgt <- mean(mean_diffs_university>mean_diff_university)
p = 2*min(c(1-pgt,pgt))
legend("topleft",legend=paste("p = ",p,sep=""),bty="n")
dev.off()

mean_diffs_DeptU <- mean_diffs_university-mean_diff_university

### Regression Analysis adjusting for Average predicted Degree
simBCArtMat <- do.call("rbind",simBCArt)
musimBCArt <- apply(simBCArtMat,2,mean)
muBCArt <- musimBCArt[which(vcols=="white")]
coefsDeptU <- numeric(1000)
coefsDeptF <- numeric(1000)
coefsNGF <- numeric(1000)
coefsUnivU <- numeric(1000)
for(i in 1:1000){
	bci <- bcArt-simBCArt[[i]][which(vcols=="white")]
	est0 <- lm(bci~ngovArt+univArt)
	coefsDeptU[i] <- coef(est0)[2]
	coefsDeptF[i] <- coef(est0)[3]
	coefsNGF[i] <- coef(est0)[2]
	coefsUnivU[i] <- coef(est0)[3]
}

pdf("BoxPlotDeptU.pdf",height=3.5,width=4,pointsize=11)
par(las=1,mar=c(3,4,1,1))
boxplot(cbind(diff = -mean_diffs_DeptU,reg = coefsDeptU),ylab="Difference in Mean Shortest Paths",ylim=c(min(c(-mean_diffs_DeptU,coefsDeptU)),max(c(-mean_diffs_DeptU,coefsDeptU))*1.25))
abline(h=0,lwd=1.5,col="red")
pgt <- c(mean(-mean_diffs_DeptU>0),mean(coefsDeptU>0))
p <- ifelse(pgt>.5,1-pgt,pgt)
legend("top",legend=paste(c("diff p = ","reg p ="),round(p,dig=4),sep=""),bty="n",horiz=T)
dev.off()

pdf("BoxPlotDeptF.pdf",height=3.5,width=4,pointsize=11)
par(las=1,mar=c(3,4,1,1))
boxplot(cbind(diff = -mean_diffs_DeptF,reg = coefsDeptF),ylab="Difference in Mean Shortest Paths",ylim=c(min(c(-mean_diffs_DeptF,coefsDeptF)),max(c(-mean_diffs_DeptF,coefsDeptF))*1.5))
abline(h=0,lwd=1.5,col="red")
pgt <- c(mean(-mean_diffs_DeptF>0),mean(coefsDeptF>0))
p <- ifelse(pgt>.5,1-pgt,pgt)
legend("top",legend=paste(c("diff p = ","reg p ="),round(p,dig=4),sep=""),bty="n",horiz=T)
dev.off()

pdf("BoxPlotNGF.pdf",height=3.5,width=4,pointsize=11)
par(las=1,mar=c(3,4,1,1))
boxplot(cbind(diff = -mean_diffs_NGF,reg = coefsNGF),ylab="Difference in Mean Shortest Paths",ylim=c(min(c(-mean_diffs_NGF,coefsNGF)),max(c(-mean_diffs_NGF,coefsNGF))*1.15))
abline(h=0,lwd=1.5,col="red")
pgt <- c(mean(-mean_diffs_NGF>0),mean(coefsNGF>0))
p <- ifelse(pgt>.5,1-pgt,pgt)
legend("top",legend=paste(c("diff p = ","reg p ="),round(p,dig=4),sep=""),bty="n",horiz=T)
dev.off()

pdf("BoxPlotUnivU.pdf",height=3.5,width=4,pointsize=11)
par(las=1,mar=c(3,4,1,1))
boxplot(cbind(diff = -mean_diffs_UnivU,reg = coefsUnivU),ylab="Difference in Mean Shortest Paths",ylim=c(min(c(-mean_diffs_UnivU,coefsUnivU)),max(c(-mean_diffs_UnivU,coefsUnivU))*2))
abline(h=0,lwd=1.5,col="red")
pgt <- c(mean(-mean_diffs_UnivU>0),mean(coefsUnivU>0))
p <- ifelse(pgt>.5,1-pgt,pgt)
legend("top",legend=paste(c("diff p = ","reg p ="),round(p,dig=4),sep=""),bty="n",horiz=T)
dev.off()


par(mfrow=c(3,3))
for(i in 1:9){
	gi <- rgraph(5)
	gplot(gi,gmode="graph",cex=4)
}










